
%% Genenerating a Worm Like Chain that is under an applied stretching force

%This code requires:
%GenTestVector.m

A=50        %nm.    persistance length
l=.1        %nm.    segment length
KT=4.2      %pN*nm
L=1000      %Bases
b=.33       %nm     length of base
N=L*b/l     %       number of segments
thetamax=pi/16;     %max angle of test vector from previous vector
showtest=1;         %show some test plots

%The main probability density function
P=inline('exp( f*l/KT*dot(t2,[0 0 1])-A/l*(1-dot(t1,t2)))','KT','f','l','A','t1','t2')

%% Testing
if showtest
    %test that I'm generating random vectors correctly
    t2=[1 1 1]'/sqrt(3);
    t=GenTestVector(t2,thetamax,500);
    tt=GenTestVector(t2,pi,500);
    
    figure;
    plot3(t(1,:),t(2,:),t(3,:),'r.')
    hold on
    plot3(tt(1,:),tt(2,:),tt(3,:),'g.')
    plot3( [0 t2(1)],[0 t2(2)],[0 t2(3)],'b','linewidth',3)
    axis equal
    axis([-1 1 -1 1 -1 1])
    
    %test that my thetamax is reasonable
    ntest=1000;
    theta=(0:ntest-1)/ntest*pi;
    for tt=1:ntest
        ttest=[sin(theta(tt)) 0 cos(theta(tt))]';
        Pcheck(tt)=P(KT,0,l,A,[0 0 1],ttest);
    end
    figure
    plot(theta,Pcheck);
    hold on
    plot([thetamax thetamax],[0 1],'r--')
    xlabel('\theta (radians)')
    ylabel('P(\theta)')
end

%% Generate Chains
forces=[0 1 10 100];
WLCs=zeros(3,N,length(forces))
for ff=1:length(forces)
    f=forces(ff);
    %build the chain using the accept-reject method
    nn=1;
    T=zeros(3,N);   %vector directions
    T(:,1)=[0 0 1]; %arbitrary starting guess
    Pmax=P(KT,f,l,0,[0 0 1],[0 0 1]); %gives max value of P
    while nn<N
        t1=T(:,nn);                     %past direction    
        t2=GenTestVector(t1,thetamax);  %test vector
        Ptest=P(KT,f,l,A,t1,t2);
        Prand=rand*Pmax;
        if Prand<Ptest %then accept
            nn=nn+1;
            T(:,nn)=t2;
            nn
        end   
    end
    
    %build chain 
    WLCs(:,:,ff)=cumsum(T,2)*l;
    
    figure
    X=WLCs(1,:,ff);
    Y=WLCs(2,:,ff);
    Z=WLCs(3,:,ff);
    plot3(X,Y,Z,'linewidth',2)
    axis equal
    xlabel('X (nm)')
    ylabel('Y (nm)')
    zlabel('Z (nm)')
    s=sprintf('Force = %g pN',f);
    title(s);
end




